Power spectrum and critical exponents in the 2D stochastic Wilson–Cowan model

The power spectrum of brain activity is composed by peaks at characteristic frequencies superimposed to a background that decays as a power law of the frequency, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f^{-\beta }$$\end{document}f-β, with an exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}β close to 1 (pink noise). This exponent is predicted to be connected with the exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}γ related to the scaling of the average size with the duration of avalanches of activity. “Mean field” models of neural dynamics predict exponents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}β and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}γ equal or near 2 at criticality (brown noise), including the simple branching model and the fully-connected stochastic Wilson–Cowan model. We here show that a 2D version of the stochastic Wilson–Cowan model, where neuron connections decay exponentially with the distance, is characterized by exponents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}β and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}γ markedly different from those of mean field, respectively around 1 and 1.3. The exponents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document}τ of avalanche size and duration distributions, equal to 1.5 and 2 in mean field, decrease respectively to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.29\pm 0.01$$\end{document}1.29±0.01 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.37\pm 0.01$$\end{document}1.37±0.01. This seems to suggest the possibility of a different universality class for the model in finite dimension.

In the last two decades, the hypothesis that the brain operates near a critical point has gained a large evidence. The first experiments pointing in this direction were done on organotypic cultures and acute slices of rat cortex 1 , where scale-free distributions of activity avalanches were found. Since then, the hypothesis has been confirmed in many systems in vitro and in vivo, from cortical activity of awake monkeys 2 to the resting MEG of human brain 3 . The distribution of avalanche sizes is found to scale as P(S) ∼ S −α , that of avalanche durations as P(T) ∼ T −τ , while the mean size of avalanches scales as �S� ∼ T γ as a function of the duration [4][5][6][7][8] . A good indicator of criticality is believed to be given by the scaling relation γ = τ −1 α−1 , as originally predicted in the theory of crackling noise 9,10 , as well as by the collapse of rescaled shapes of avalanches of different durations 4,7,8 . The simple branching model of avalanche propagation predicts exponents α = 3/2 , τ = 2 and γ = 2 , observed in some experimental realizations 1,11 and in models of neural dynamics, including the fully-connected stochastic Wilson-Cowan model 7 . However, some experimental results have found an exponent γ around 1.3, not compatible with the value 2 predicted by the branching process universality class, even when the scaling relation is experimentally satisfied. For instance, results on spike avalanches measured in the urethane-anesthetized rat cortex 12 , cultured cortical networks 4 , ex-vivo recordings of the turtle visual cortex 5 , and somatosensory barrel cortex of the anesthetized rat 13 , have found γ around 1.3 with the scaling relation γ = τ −1 α−1 satisfied. Another important feature of neuronal dynamics is the power-law decay P(f ) ∼ f −β of the power spectrum of EEG, MEG, resting state fMRI, and local field potential as a function of frequency [14][15][16][17][18] , once the peaks corresponding to characteristic frequencies of oscillations have been subtracted. The values of β are found to be between 1 and 1.3 in the EEG and MEG of healthy patients 16,17 , while they are around 2 for epileptic patients 19 . On quite general grounds, the exponent β is predicted to be equal to the exponent γ of the relation �S� ∼ T γ9 . It seems therefore that there are at least two universality classes in brain dynamics, one that can be called the "mean-field" class, represented by the branching model and the fully-connected Wilson-Cowan model, that is characterized by α ≃ 3/2 and β ≃ γ ≃ τ ≃ 2 (the equality is exact for the branching model), and another characterized by a lower value of the exponents, β ≃ γ 1.3 , α ≃ 1.3 , τ ≃ 1.4 . The experimental measurement of γ around 1.3, not compatible with the value of branching process, opened a debate, and it raised the question if a different model, with a different universality class 12,20 , might be required to explain the critical brain data 4,5,12,13 , or if different mechanisms, such as subsampling 21 , should be invoked to reconcile data with branching model. www.nature.com/scientificreports/ In the present paper we study the 2D version of the stochastic Wilson-Cowan model 7,22 , where neurons are distributed uniformly on a 2D lattice, and connections between them decay exponentially with the distance. We can tune the network topology from 2D to fully-connected by changing the ratio between the range of the exponential decay of the connections, and the side of the lattice L. While in the fully-connected case the dynamics of the model is completely described by specifying two variables, the fraction of active excitatory and inhibitory neurons, in the 2D case one needs to define such fractions at each site of the lattice. As a consequence, while the fully-connected model is characterized by just two characteristic relaxation times, in the 2D model one finds a whole spectrum of times, related to the different Fourier modes of the neural activity. By making a system size expansion one finds that, if the number of neurons on each site of the lattice is large, the dynamical equations governing the evolutions of the Fourier modes decouple, and each mode obeys to the same equations of the fully-connected model, but with a different relaxation time.
Independently of the network topology, the model shows a critical point at a characteristic value of a parameter measuring the difference between the strength of excitatory and inhibitory connections. For values above the critical point, the model displays a self-sustained dynamics even in the absence of external input. At the critical point, one of the characteristic times diverges, related to the mode k = 0 in spatially extended topologies, while the times of the other modes scale as |k| −2 . This feature, together with the density of the wave numbers which scales as |k| in two dimensions, gives rise to a power spectrum that is proportional to f −1 (pink noise).
In the next section we introduce the model, then we study the power spectrum and relaxation functions of the firing rate in the linear case (large neuron density). We then look at the avalanche size and duration distributions, and show that at the critical point the distributions are scale free, with exponents depending on the topology of the system. Finally, we show that the exponents β and γ have approximately the same dependence on the inverse frequency and duration of avalanches, respectively, and in the 2D model tend respectively to 1 and 1.3 at low frequency and large durations.

The model
Let us consider a two-dimensional L × L lattice, where on each site there are n E = N E /L 2 excitatory and n I = N I /L 2 inhibitory neurons, with connections depending on the distance r ij measured in lattice spacings, as shown in the scheme in Fig. 1. Note that all pairs i and j of neurons belonging to the same site have distance r ij = 0 , so that inside one lattice site the network is fully-connected. Neurons are modeled as in 22 , namely they can be in two states, active and quiescent. The rate of transition from active to quiescent state is α for all the neurons, while the rate from quiescent to active state is given by an activation function f (s i ) of the input s i of the i-th neuron, given by where a j = 0, 1 if the j-th neuron is quiescent or active respectively, w ij are the connections between neurons, and h i is an external input. We consider the activation function In the following we set α = 0.1 ms −1 , and β = 1 ms −1 . We study here a version of the model where the connections between neurons do not depend only on the type of presynaptic neuron, as in the fully-connected case, but depend also on the distance between neurons. Namely, the connection between neurons j (pre-synaptic) and i (post-synaptic) is given by (3) www.nature.com/scientificreports/ where we have defined w 0,0 = j w ij as the sum of all the connections incoming in one neuron, which we take to be the same for all neurons, w s,0 = j |w ij | , and r ij is the distance between neurons. The second subscript in w 0,0 and w s,0 refers to the wave-number k = 0 , see Eq. (13). The normalization factors N E and N I are defined as N E = j∈E e −r ij / , N I = j∈I e −r ij / , where E and I are respectively the set of all excitatory and inhibitory neurons. Note that, when ≫ L , we recover the fully-connected model, well studied in the past 7,22,23 . On the other hand, when L ≫ , the topology of the connections changes to two-dimensional. We consider the same input h i ≡ h for all neurons. As the connections w ij depend only on the distance between neurons, and the input h is the same for different neurons, the system is translationally invariant. The fraction of active neurons is a stochastic variable that at stationarity fluctuates around the fixed point value, given by the equation (see "Methods" section) with s 0 =w 0,0 0 + h . Note that Eq. (4) does not depend on the chosen topology of the connections, fully-connected or sparse or depending on the distance, but only on the condition that the sum of incoming connections w 0,0 is the same for all neurons, because this makes the system translationally invariant and the fixed point activity 0 equal for all the lattice sites. In Ref. 7 it was shown that there is a critical point at h = 0 and w 0,0 = β −1 α . For w 0,0 larger than the critical value, an attractive fixed point with � 0 > 0 exists even when the external input h → 0 . In the case of the fully-connected model, the connections w ij do not depend on the spatial position of neurons, but only on the functional type (excitatory or inhibitory) of the pre-synaptic and post-synaptic neuron. In Refs. 7,22,23 a further simplification was considered, that w ij depends only on the type of the pre-synaptic neuron, and is given by In this case, the temporal autocorrelation function of time dependent variables, like the fraction of active neurons or the firing rate, can be written in the limit of large number of neurons as 23 where X(t) is the variable considered, X 0 is its average value in time, A 1 and A 2 are constant coefficients, and τ 1 , τ 2 are the characteristic relaxation times. While τ 2 is always small, and lower than α −1 , the time τ 1 diverges at the critical point 7 .
In the fully-connected case, the state of the system depends only on two variables and , that correspond to the sum and difference between the fractions of activated excitatory and inhibitory neurons; conversely for connections depending on the distance, the values r and r on each site of the lattice are necessary to characterize the activity. They are defined as where m r and l r are respectively the number of active excitatory and inhibitory neurons on the lattice site r . Equivalently, we can use the Fourier transforms k and k , where k 's are L 2 different wave vectors. As the system is translationally invariant, the fixed point is characterized by k = k = 0 for k = 0 . Moreover, for the wave vector k = 0 , at the fixed point 0 = 0 (this is a consequence of the fact that connections do not depend on the post-synaptic neuron), while 0 obeys the same Eq. (4) of the fully-connected case. Therefore, for w 0,0 ≤ β −1 α , the fixed point value 0 goes to zero when the external input h → 0 , while for w 0,0 > β −1 α it remains finite even for h → 0.
A quantity that will be considered in the following is the local firing rate, defined as where s r is the input defined by Eq. (1) of neurons on site r . Note that all the neurons on site r have the same input. The probability that a neuron on site r fires (makes a transition from inactive to active state) in the interval of time t is R r (t)�t.

Temporal correlations and power spectrum
We now consider a variable X r (t) defined on a site r , that can be for instance the fraction of active neurons � r (t) , or the firing rate R r (t) , given by a superposition of all its Fourier components, As shown in "Methods", for large number of neurons different Fourier components of the activity decouple, and obey the same evolution equation of the total activity in the system with full connectivity. Therefore the autocorrelation function of X r (t) will be given by the sum of the autocorrelation functions of its Fourier components, www.nature.com/scientificreports/ While A 2,k and τ 2,k remain always finite and small, A 1,k and τ 1,k diverge proportionally to |k| −2 at the critical point. (17). The power spectrum of the activity on a site of the lattice is given by the Wiener-Khinchin theorem as the temporal Fourier transform of the autocorrelation function, and given the dependence of A 1,k and τ 1,k on the wave number at the critical point, in a system with a two-dimensional structure one finds a dependence P(f ) ∝ 1/f for frequencies between f min = L −2 D 0 and f max = −2 D 0 , where L is the linear size of the lattice and the range of the connections (see "Methods"). Note that the exponent of the power spectrum depends on the dimension of the space, in particular . At frequencies f ≫ f max only the Fourier components with the fastest relaxation time f −1 max survive in the spectrum, so that P(f ) ∝ f −2 in any spatial dimension.
In Fig. 2A,B we show the dependence of the autocorrelation function and of the power spectrum on the distance with respect to the critical point, which is given by w 0,0 = β −1 α = 0.1 and h = 0 . Far from the critical point, the distribution of times for different wave vectors is narrow, the autocorrelation decays exponentially with good approximation and the power spectrum can be described by a Lorentzian. Near the critical point, we have a wide distribution of times, that gives rise to a 1/f decay in the spectrum. As shown in Fig. 2A, the corresponding relaxation function exhibits the slow decays a − b ln t for a wide interval of times 24 .
As shown in Fig. 3A,B, the 1/f dependence of the spectrum depends not only on the distance from the critical point, but also on the range of the connections. If the range grows, the autocorrelation functions and power spectrum tend to those of the fully-connected system, that is characterized by just two correlation times τ 1 and τ 2 , where τ 2 is a short time of the order of α −1 while τ 1 is large near the critical point, so that the autocorrelation is well described at long times by a single exponential. In Fig. 3C, we show the power spectrum evaluated in systems where connections decay as a power law. Namely, in Eq. (3), instead of the factor e −r/ we put a factor min(1, r −ω ) . The case ω = 0 coincides with the case = ∞ , and gives again the 1/f 2 decay of the spectrum. Larger values of ω correspond to connections decaying more quickly, and for ω ≥ 4 one finds the 1/f decay characteristic of short range connections.

Relation between avalanche and power spectrum exponents
In this section we compare the size and duration distributions of avalanches of activity of a site of the lattice in the 2D case (with short range connections) and in the fully-connected model. In the following we set n E = n I = 10 8 neurons for each site, and compare the behaviour of the network with L = 40 and = 1 (2D), with that of a network with L = 1 (fully-connected). We remark here that in a network with L = 1 the distance between all neurons is zero, so that in Eq. (3) the connections w ij depend only on the type of neuron, and the model coincides with the model studied in 7,22,23 , with full connectivity. We simulated the system with Langevin dynamics, see Eq. (7). To speed up the simulations, we have set the connections w ij in Eq. (3) to zero when r ij > 3 . We made 60 different runs both for L = 1 and L = 40 , for respectively 3.5 × 10 7 ms and 2.5 × 10 5 ms, discarding the first 10 7 and 4 × 10 4 ms, and collecting the avalanches on all sites in the 2D system. We define an avalanche as follow: we divide the time in discrete bins of width δ = 1 ms, and consider the time evolution of the activity on a single site. We identify an avalanche as a continuous series of time bins in which there is at least one spike (i.e., a transition of one neuron from a quiescent to an active state www.nature.com/scientificreports/ in the site being considered). Note that, when simulating the system with Langevin dynamics Eq. (7), we extract the number of spikes in the interval δ from a Poissonian distribution with mean equal to the temporal integral over δ of the firing rate defined in Eq. (5) multiplied by the number of neurons in the site. We remark here that avalanches are relative to the activity on a single site, so an avalanche begins and ends when the activity on the considered lattice site is zero, regardless of the activity on the other sites of the lattice. The size of the avalanche is defined as the total number of spikes of neurons belonging to the site considered, while the duration is the number of time bins of the avalanche multiplied by the width δ of the bins. Note that all the lattice sites are equivalent, because the connections are translationally invariant and boundary conditions are periodic, so we expect the same distribution of sizes and durations on all the sites of the network. Therefore, to improve statistics, we compute the average of the distributions over all the sites of the lattice.
The distributions of avalanche size and duration are expected to follow power-law scalings, P(S) ∼ S −α and P(T) ∼ T −τ . In Fig. 4 we report the results for the 2D system ( L = 40 , blue dots) and the fully-connected system ( L = 1 , red dots), for the size distribution (A) and duration distribution (B). The distributions show a clear dependence on the topology of the network. In the fully-connected model (red dots), exponents are α = 1.48 ± 0.01 and τ = 2.05 ± 0.01 , very close to the values characteristic of the branching model of neural dynamics, as already observed in 7 . In the 2D system (blue dots), after a small size and duration regime characterized by exponents close to one, one finds, up to the cut-off of the distributions, a scaling regime with exponents α = 1.29 ± 0.01 and τ = 1.37 ± 0.01 , markedly smaller than the mean-field values.
We have then considered the relation S (T) between the duration of an avalanche and its mean size (Fig. 5A), for the same system and parameters of Fig. 4. We fit the measured data, in the same ranges of Fig. 4B, with a power law T γ . Also in this case, we observe a marked difference between the fully-connected model (red dots), where the exponent is γ = 2.04 ± 0.01 , and the 2D system where γ = 1.35 ± 0.01 . Note that in both cases the expected relation 9,10   To investigate the relation between the exponent γ of S (T) and the exponent β of the power spectrum, in Fig. 5D we plot the power spectrum P(f) of the single site firing rate defined in Eq. (5), for the same system size and parameters of Fig. 5A,B and Fig. 4. The power spectrum of the fully-connected system is well fitted by a Lorentzian, with a white noise behaviour for frequencies lower than 1 Hz, and a decay with an exponent β = 1.98 ± 0.02 for frequencies larger than 1 Hz. Conversely the spectrum of the 2D system, as anticipated analitically for the N → ∞ case, shows a decay with an exponent β = 1.02 ± 0.02 for frequencies between 0.05 and 1 Hz, intermediate between white noise at lower frequencies (not shown) and brown noise at higher ones. As we have done with the exponent γ in Fig. 5B, in Fig. 5E we plot the exponent β of a fit P(f ) ∼ f −β in a sliding window of range [0.1 f , f ] . Comparing Fig. 5E with Fig. 5B, it can be seen that exponents γ and β have a similar dependence respectively on the avalanche duration and on the inverse frequency. In the case of the fully-connected system (red dots) both exponents are around 2 in the ranges [10 2 , 10 3 ] ms and [1,10] Hz, and decay to lower values for longer avalanches or smaller frequencies. In the case of the 2D system the exponents are nearly constant in the ranges [10 3 , 10 4 ] ms and [0.1, 1] Hz, although they tend to quite different values, namely β ≃ 1 , γ ≃ 1.3 . The reason of this discrepancy is to be further investigated.
In Fig. 5C we show S (T) for different sizes of the lattice, from L = 1 (fully-connected case) to L = 40 , while in Fig. 5F we plot the exponent γ measured by fits of the data in panel C in a range of mean sizes between 5 × 10 6 and 10 8 . It can be seen that the exponent remains near to γ = 2 for small sizes, and drops to a value near γ = 1.3 for large lattices.

Scaling of the shape of avalanches
Together with the relation Eq. (6), another test of the "criticality" hypothesis for avalanche activity, is the scaling of the shape of the avalanches. Denoting with V(t) the mean number of spikes observed at time t during an avalanche of duration T, the total size of the avalanche is given by S = T 0 dtV (t) ∼ T γ , so it is expected that the normalized shape V (t)T 1−γ depends only on the rescaled time t/T. This relation should hold as long as the mean size S (T) is well fitted by a power law T γ , that is for durations where the exponent of the sliding window fits are nearly constant in Fig. 5B. Looking at Fig. 5B, we expect a collapse of the shapes in the interval 10 2 < T < 10 3 ms for the fully-connected system, with a value γ ∼ 2.04 , and in the interval 10 3 < T < 10 4 ms for the 2D system, with a value γ ∼ 1.35 . We highlight such values of γ by dashed lines in the figure. The collapse of the shapes is indeed quite well verified, for such values and duration ranges, as shown in Fig. 6A,B. Note that, as expected, an exponent γ ≃ 2 corresponds to a shape that is nearly parabolic (fully-connected system), while in the case of the 2D system the exponent γ ≃ 1.3 corresponds to a more flattened shape 26 .

Dependence on the fraction of inhibitory neurons
In this section we investigate the role of the ratio between excitatory and inhibitory neurons, by simulating the system also for a fraction 80/20 of excitatory and inhibitory neurons, more similar to cortical networks. The total number of neurons on each site of the lattice is 2 × 10 8 as in the previously studied (50/50) case, so that we have n E = 1.6 × 10 8 and n I = 4 × 10 7 . Note that, as it is apparent from Eq. (3), the strength of the excitatory/inhibitory synapses is inversely proportional to the number of excitatory/inhibitory neurons, so that if the number of inhibitory neurons is decreased, the strength of inhibitory connections is correspondingly increased. In this way we ensure that the critical point corresponds to the same value of w 0,0 = β −1 α and h = 0.
In Fig. 7 we compare the four cases considered: fully-connected 50/50, fully-connected 80/20, 2D 50/50 and 2D 80/20. We note that the fraction of inhibitory neurons affects only the cut-off in the distributions P(S) and www.nature.com/scientificreports/ P(T), while the mean avalanche size S (T) and the power spectrum P(f) are not affected. Moreover the cut-off has a significative change only in the 2D system, and is much less affected in the fully-connected case. In the insets of Fig. 7A,B, we report a collapse of the 50/50 and 80/20 distributions in the 2D case, showing that the cutoff is approximately four (three) times smaller in the 80/20 case for the size (duration) distributions.
In conclusion, we can affirm that the observed dependence of critical exponents on the spatial dimensionality is preserved also for different fractions of inhibitory neurons.

Conclusions
We have studied the stochastic Wilson-Cowan model on a 2D lattice, with connections decaying exponentially with the distance. We use a connection weight that decays exponentially with distance to model the structural anatomical connectivity that exhibits exponential decay with distance. Recent research indeed, using data from retrograde tracer injections, shows that influence of interareal distance on connectivity patterns is conform to an exponential distance rule (EDR), according to which the projection lengths decay exponentially with interareal separation [28][29][30][31] . Spatial decay constant is of order of few mm for white matter, and few tenths of mm for gray matter, in mouse and macaque, 28,29 . Notably, in the 2D model, as in the fully-connected case, varying the difference between the strength of excitatory and inhibitory connections, the model undergoes a second order transition from a phase where the activity tends to zero in absence of external input, to a phase where the activity is self-sustained even in absence of external inputs. At the critical point, one of the relaxation times diverge. While the fully-connected model is characterized by just two relaxation times, one of which is always lower than α −1 , where α is the deactivation rate of neurons, the 2D system has a spectrum of relaxation times, related to the different Fourier modes that can be defined on the lattice. Such relaxation times become proportional to D −1 0 |k| −2 at the critical point, where D 0 is a constant and, as a consequence, the model in 2D shows a logarithmic decay of the relaxation functions and a f −1 behaviour of the power spectrum, in an interval of frequencies between f min = L −2 D 0 and f max = a −2 D 0 , where L is the lattice size, is the range of connections, and D 0 is defined in Eq. (17). We emphasize that the 1/f behaviour is observed for the single site activity (or the activity of a localized group of sites), being the superposition of all the Fourier components with a spectrum of relaxation times. Conversely, the activity of the entire system corresponds to the single Fourier component k = 0 and therefore exhibits www.nature.com/scientificreports/ a 1/f 2 behavior. Moreover, we have shown that the change in the behaviour of the power spectrum is related to a marked change in the exponents of the avalanche distributions. Note that, in the case of the 2D system, we have defined avalanches considering the activity on a single site and not in the whole system, in order to use the same signal considered to define the power spectrum. This choice has been also determined by the observation that, in a large 2D system, activity in the entire system never goes to zero, which would make the introduction of a thresholding procedure necessary to define avalanches. Clearly, although avalanches are measured locally, their distribution reflects the activity of all the sites of the system. In the 2D case, the exponents α and τ of the size and duration distributions become smaller than the ones predicted in the mean-field case, respectively α ≃ 1.3 and τ ≃ 1.4 , while the exponent γ that relates the mean size to the avalanche duration decreases from 2 to 1.3. It remains to be understood the discrepancy between γ and the exponent β of the power spectrum. Notably in the 2D model the avalanche shape collapse is quite well verified with a value γ ≃ 1.35 , showing a more flattered shape, while it is γ ≃ 2 in the fully-connected case. Discrepancy between the fully-connected γ ≃ 2 prediction and the exponent γ ≃ 1.3 observed in some experiments 12 , with a large interval of exponents α and τ of the size and duration distributions all falling on the γ ≃ 1.3 line, may be attributed to the effective topology of the measured activity. For example measurement able to record spiking activity of large areas of neuronal populations may reflect the structured topology of the network, while more localized measures (in highly connected areas) may be better approximated by mean-field i.e. fully-connected topology. It would be interesting to compare the corresponding scaling behavior of the power spectra of experimental activity in different topological conditions, to confirm the existence of different universality classes.

Methods
We consider a two-dimensional lattice of L × L sites. On each site there are n E = N E /L 2 excitatory neurons and n I = N I /L 2 inhibitory ones. Define m r and l r respectively the number of active excitatory and inhibitory neurons on the site r = (x, y) , with x, y = 0, . . . , L − 1.
In the Gaussian noise approximation, the temporal evolution of the system can be effectively described in terms of the coupled non-linear Langevin equations 32 where α is the rate of deactivation of the neurons, f(s) is the input dependent rate of activation, s  www.nature.com/scientificreports/ We can perform a system size expansion, changing variables from the fractions of active neurons x r and y r to the deviations of the fractions with respect to the fixed point. Defining x r = x 0 + n −1/2 E ξ E,r , y r = y 0 + n −1/2 I ξ I,r , substituting in Eq. (8), and neglecting terms that are small when the number of neurons is large, we obtain the linear Langevin equations The Eq. (9) are translationally invariant and linear, therefore it is convenient to perform a Fourier transform and define where k = 2π L (k x , k y ) , with k x , k y = 0, . . . , L − 1. Substituting in (9), we obtain the decoupled equations for each of the L × L pairs of Fourier modes ξ E,k and If the connections do not depend on the post-synaptic neuron, but only on the pre-synaptic one, that is w r ′ r . Note that τ 2,k is independent of k . The fixed point input s 0 can be written as s 0 =w 0,0 0 + h , and we report here for convenience the form of the fixed point equation, As we have seen, each of the Fourier modes behaves exactly as the total fraction of neurons in the model with all-to-all connections, but with parameters τ 1,k , w ff,k that depend on the Fourier mode k . The autocorrelation function of the fluctuations ξ �,k is therefore given by 23 rr ′ ξ I,r ′ + 2αy 0 η I,r (t).
�ξ �,k (t)ξ �,−k (0)� = A 1,k e −t/τ 1,k + A 2,k e −t/τ 2,k , www.nature.com/scientificreports/ By performing an inverse Fourier transform, we can find the autocorrelation function of ξ �,r = ξ E,r +ξ I,r 2 , which is given by The power spectrum of the variable ξ �,r (t) is therefore given by the Wiener-Khinchin theorem as The linear Eq. (10) depend basically on parameters defined at the fixed point. The fixed point can undergo different kinds of transitions (bifurcations). We will consider here the case of the transcritical bifurcation, where one of the eigenvalues of the matrix in Eq. (11) vanishes, in particular we will consider the case in which the eigenvalue τ −1 1,0 corresponding to the mode k = 0 vanishes, so that For values of |k| smaller than −1 , where is the range of the connections w The values of τ 1,k therefore diverge for |k| → 0 , so that we can consider τ 1,k ≫ τ 2,k for low values of the wave number. In this case, we can approximate Therefore both τ 1,k and A 1,k diverge for |k| → 0 , while τ 2,k , A 2,k e w ff,k remain finite. Neglecting terms relative to τ 2,k in the power spectrum, we obtain In d spatial dimensions, the density of wave numbers is L Note that the previous derivation holds only for ≪ L , otherwise the condition f min ≪ f ≪ f max cannot be met.

Data availability
The datasets analysed during the current study are available from the corresponding author on reasonable request.